second approach

#library(R2jags)
library(R2jags)
library(tidyverse)
library(lubridate)
library(knitr)
library(broom)

leaders_data <- get(load("Leaders.RData"))

leaders_data$Censored <- as.factor(leaders_data$Censored)
leaders_data$Type <- as.factor(leaders_data$Type)
leaders_data$Age.Event <- as.numeric(leaders_data$Age.Event)
summary(leaders_data)

##      Name             Birth.Date           Event.Date           Age.Event     
##  Length:177         Min.   :1324-05-14   Min.   :1368-03-29   Min.   : 9.259  
##  Class :character   1st Qu.:1507-09-16   1st Qu.:1572-07-05   1st Qu.:53.985  
##  Mode  :character   Median :1705-10-31   Median :1758-05-03   Median :67.529  
##                     Mean   :1674-01-23   Mean   :1737-03-27   Mean   :63.170  
##                     3rd Qu.:1833-08-20   3rd Qu.:1886-11-18   3rd Qu.:78.242  
##                     Max.   :1961-08-04   Max.   :2020-07-31   Max.   :95.830  
##  Censored        Type         Fail           centdate               time       
##  0:166    ChinaEmp :26   Min.   :0.0000   Min.   :1300-01-01   Min.   :0.2436  
##  1: 11    DalaiLama:14   1st Qu.:1.0000   1st Qu.:1300-01-01   1st Qu.:2.0770  
##           JapanEmp :30   Median :1.0000   Median :1300-01-01   Median :4.0582  
##           Pope     :63   Mean   :0.9379   Mean   :1300-01-01   Mean   :3.7406  
##           UsPres   :44   3rd Qu.:1.0000   3rd Qu.:1300-01-01   3rd Qu.:5.3362  
##                          Max.   :1.0000   Max.   :1300-01-01   Max.   :6.6157

Boxplot of age (Age.Event) by Type

ggplot(data = leaders_data) +
  geom_boxplot(aes(x = Type, y = Age.Event))

ggplot(data = leaders_data) +
  geom_point(aes(x = Birth.Date, y = Age.Event))

THIS IS OUR MODEL

\[ T_i \sim Weibull(\alpha, \mu) \\ log(\mu) = \beta_0 + \beta_1x_{i1}+\beta_2x_{i2}+...+\beta_px_{ip}\\ \] Log of mu is done so that the mu parameter is always positive (the support on mu in Weibull is that mu > 0)

More specificically

\[ T_i \sim Weibull(r, \mu_i) \\ \log(\mu_i) = \beta_0+\beta_1 \text{(Year of Birth_i)}+\beta_2 \text{(Leadership_i = UsPres)} + \beta_3\text{(Leadership_i = ChinaEmp)} \\ + \beta_4\text{(Leadership_i = DalaiLama)} +\beta_5\text{(Leadership_i = JapanEmp)} \\ + \beta_6\text{(Year of Birth_i * (Leadership_i = UsPres))} \\ + \beta_7\text{(Year of Birth_i * (Leadership_i = ChinaEmp))} \\ + \beta_8\text{(Year of Birth_i * (Leadership_i = DalaiLama))} \\ + \beta_9\text{(Year of Birth_i * (Leadership_i = JapanEmp))} \]

Data Cleaning

leaders_data$Birth.Year <- year(leaders_data$Birth.Date)

leaders_data <- leaders_data %>%
  mutate(TypeChinaEmp = as.factor(case_when(Type == "ChinaEmp" ~ 1,
                             TRUE ~ 0)),
         TypeDalaiLama = as.factor(case_when(Type == "DalaiLama" ~ 1,
                                   TRUE ~ 0)),
         TypeJapanEmp = as.factor(case_when(Type == "JapanEmp" ~ 1, 
                                  TRUE ~ 0)),
         TypePope = as.factor(case_when(Type == "Pope" ~ 1, 
                              TRUE ~ 0)),
         TypeUsPres = as.factor(case_when(Type == "UsPres" ~ 1,
                                TRUE ~ 0)))

leaders_nopred <- leaders_data %>%
  filter(!(Name %in% c("14th Dalai Lama (Tenzin Gyatso)", "Naruhito (Emperor Reiwa)", "Barack Obama", "Francis")))

survival_raw <- leaders_nopred$Age.Event
n <- length(survival_raw)

censored <- leaders_nopred$Censored
survival <- ifelse(censored == 0, survival_raw, NA)

censoring_limits <- ifelse(censored == 0, 120, survival_raw)
x_1 <- leaders_nopred$Birth.Year - mean(leaders_nopred$Birth.Year)
x_2 <- leaders_nopred$TypeUsPres
x_3 <- leaders_nopred$TypeChinaEmp
x_4 <- leaders_nopred$TypeDalaiLama
x_5 <- leaders_nopred$TypeJapanEmp
# add interaction b/w Birth.Year and Type
x_6 <- x_1*as.numeric(as.character(x_2))
x_7 <- x_1*as.numeric(as.character(x_3))
x_8 <- x_1*as.numeric(as.character(x_4))
x_9 <- x_1*as.numeric(as.character(x_5))

set.seed(123)
second_approach <- function() {
  for(i in 1:n_censored) {
    z_censored[i] ~ dpois(phi_censored[i])
    phi_censored[i] <- mu_censored[i] * pow(t_censored[i], r)
    mu_censored[i] <- exp(beta_censored[i])
    beta_censored[i] <- beta_0 + beta_1*x_1_censored[i] + beta_2*x_2_censored[i] + beta_3*x_3_censored[i] + beta_4*x_4_censored[i] + beta_5*x_5_censored[i] + beta_6*x_6_censored[i] + beta_7*x_7_censored[i] + beta_8*x_8_censored[i] + beta_9*x_9_censored[i]
  }
  
  for(j in 1:n_non_censored) { #total rows - 7 ** 7 are censored, rest are not censored
    survival_non_censored[j] ~ dweib(r, mu[j])
    mu[j] <- exp(beta[j])
    beta[j] <- beta_0 + beta_1*x_1_non_censored[j] + beta_2*x_2_non_censored[j] + beta_3*x_3_non_censored[j] + beta_4*x_4_non_censored[j] + beta_5*x_5_non_censored[j] + beta_6*x_6_non_censored[j] + beta_7*x_7_non_censored[j] + beta_8*x_8_non_censored[j] + beta_9*x_9_non_censored[j]
  }
  
  beta_0 ~ dnorm(0.0, 1.0E-3) # Prior on beta_0 is normal with low precision
  beta_1 ~ dnorm(0.0, 1.0E-3) # Prior on beta_1 is normal with low precision
  beta_2 ~ dnorm(0.0, 1.0E-3) # Prior on beta_2 is normal with low precision
  beta_3 ~ dnorm(0.0, 1.0E-3) # Prior on beta_3 is normal with low precision
  beta_4 ~ dnorm(0.0, 1.0E-3) # Prior on beta_4 is normal with low precision
  beta_5 ~ dnorm(0.0, 1.0E-3) # Prior on beta_5 is normal with low precision
  beta_6 ~ dnorm(0.0, 1.0E-3) # Prior on beta_6 is normal with low precision
  beta_7 ~ dnorm(0.0, 1.0E-3) # Prior on beta_7 is normal with low precision
  beta_8 ~ dnorm(0.0, 1.0E-3) # Prior on beta_8 is normal with low precision
  beta_9 ~ dnorm(0.0, 1.0E-3) # Prior on beta_9 is normal with low precision
  r ~ dexp(1) # Prior on r
  
  # Predictive distribution of age at the new values
  beta_Francis <- beta_0 + beta_1*year_Francis
  mu_Francis <- exp(beta_Francis)
  survival_Francis ~ dweib(r, mu_Francis) %_% T(present_length_Francis, upper_length)
  age_Francis_predictive <- survival_Francis

  beta_Obama <- beta_0 + (beta_1+beta_6)*year_Obama + beta_2
  mu_Obama <- exp(beta_Obama)
  survival_Obama ~ dweib(r, mu_Obama) %_% T(present_length_Obama, upper_length)
  age_Obama_predictive <- survival_Obama

  beta_Dalai <- beta_0 + (beta_1+beta_8)*year_Dalai + beta_4
  mu_Dalai <- exp(beta_Dalai)
  survival_Dalai ~ dweib(r, mu_Dalai) %_% T(present_length_Dalai, upper_length)
  age_Dalai_predictive <- survival_Dalai
  
  beta_Naruhito <- beta_0 + (beta_1+beta_9)*year_Naruhito + beta_5
  mu_Naruhito <- exp(beta_Naruhito)
  survival_Naruhito ~ dweib(r, mu_Naruhito) %_% T(present_length_Naruhito, upper_length)
  age_Naruhito_predictive <- survival_Naruhito
}

censored_dex = which(censored==1) # index of ppl who are censored

z_censored <- rep(0, length(censored_dex))
t_censored <- censoring_limits[censored_dex]
x_1_censored <- x_1[censored_dex] 
x_2_censored <- x_2[censored_dex] 
x_3_censored <- x_3[censored_dex]
x_4_censored <- x_4[censored_dex]
x_5_censored <- x_5[censored_dex]
x_6_censored <- x_6[censored_dex]
x_7_censored <- x_7[censored_dex]
x_8_censored <- x_8[censored_dex]
x_9_censored <- x_9[censored_dex]

n_censored <- length(censored_dex)

survival_non_censored <- survival[-censored_dex] # Remove ALL ALIVE PEOPLE
x_1_non_censored <- x_1[-censored_dex]
x_2_non_censored <- x_2[-censored_dex]
x_3_non_censored <- x_3[-censored_dex]
x_4_non_censored <- x_4[-censored_dex]
x_5_non_censored <- x_5[-censored_dex]
x_6_non_censored <- x_6[-censored_dex]
x_7_non_censored <- x_7[-censored_dex]
x_8_non_censored <- x_8[-censored_dex]
x_9_non_censored <- x_9[-censored_dex]

n_non_censored <- length(survival_non_censored)

birth_year_Francis <- leaders_data$Birth.Year[leaders_data$Name == "Francis"]
year_Francis <- birth_year_Francis - mean(leaders_nopred$Birth.Year)
present_length_Francis <- leaders_data$Age.Event[leaders_data$Name =="Francis"]

birth_year_Obama <- leaders_data$Birth.Year[leaders_data$Name == "Barack Obama"]
year_Obama <- birth_year_Obama - mean(leaders_nopred$Birth.Year)
present_length_Obama <- leaders_data$Age.Event[leaders_data$Name =="Barack Obama"]

birth_year_Dalai <- leaders_data$Birth.Year[leaders_data$Name == "14th Dalai Lama (Tenzin Gyatso)"]
year_Dalai <- birth_year_Dalai - mean(leaders_nopred$Birth.Year)
present_length_Dalai <- leaders_data$Age.Event[leaders_data$Name =="14th Dalai Lama (Tenzin Gyatso)"]

birth_year_Naruhito <- leaders_data$Birth.Year[leaders_data$Name == "Naruhito (Emperor Reiwa)"]
year_Naruhito <- birth_year_Naruhito - mean(leaders_nopred$Birth.Year)
present_length_Naruhito <- leaders_data$Age.Event[leaders_data$Name =="Naruhito (Emperor Reiwa)"]

upper_length <- 120

data_Popes_second_approach <- list("n_censored",
                                   "n_non_censored",
                                   "z_censored",
                                   "t_censored",
                                   "x_1_censored",
                                   "x_2_censored",
                                   "x_3_censored",
                                   "x_4_censored",
                                   "x_5_censored",
                                   "x_6_censored",
                                   "x_7_censored",
                                   "x_8_censored",
                                   "x_9_censored",
                                   "survival_non_censored",
                                   "x_1_non_censored",
                                   "x_2_non_censored",
                                   "x_3_non_censored",
                                   "x_4_non_censored",
                                   "x_5_non_censored",
                                   "x_6_non_censored",
                                   "x_7_non_censored",
                                   "x_8_non_censored",
                                   "x_9_non_censored", 
                                   "year_Francis",
                                   "year_Obama", 
                                   "year_Dalai",
                                   "year_Naruhito",
                                   "present_length_Francis", 
                                   "present_length_Obama",
                                   "present_length_Dalai",
                                   "present_length_Naruhito",
                                   "upper_length")

Bayesian_Popes_second_approach <- jags(data = data_Popes_second_approach,  
                                      parameters.to.save = c("beta_0",
                                                             "beta_1",
                                                             "beta_2", 
                                                             "beta_3", 
                                                             "beta_4", 
                                                             "beta_5",
                                                             "beta_6",
                                                             "beta_7", 
                                                             "beta_8", 
                                                             "beta_9",
                                                             "r", 
                                                             "age_Francis_predictive",
                                                             "age_Obama_predictive",
                                                             "age_Dalai_predictive",
                                                             "age_Naruhito_predictive",
                                                             "survival_non_censored"
                                                             ), 

                                      n.iter = 50000, 
                                      n.chains = 3,
                                      model.file = second_approach)

## module glm loaded

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 173
##    Unobserved stochastic nodes: 15
##    Total graph size: 2385
## 
## Initializing model

Bayesian_Popes_second_approach

## Inference for Bugs model at "/tmp/RtmpyN0MZP/model49d3f84c442.txt", fit using jags,
##  3 chains, each with 50000 iterations (first 25000 discarded), n.thin = 25
##  n.sims = 3000 iterations saved
##                             mu.vect sd.vect     2.5%      25%      50%      75%
## age_Dalai_predictive        100.942   9.940   85.743   92.172  100.329  109.281
## age_Francis_predictive      104.921  10.069   85.599   96.968  106.165  113.568
## age_Naruhito_predictive      98.891  15.455   66.047   87.555  101.833  112.140
## age_Obama_predictive         99.278  15.361   64.529   89.100  102.630  112.128
## beta_0                      -22.776   1.596  -25.941  -23.801  -22.745  -21.742
## beta_1                       -0.002   0.001   -0.003   -0.002   -0.002   -0.001
## beta_2                        0.601   0.440   -0.282    0.307    0.610    0.906
## beta_3                        1.324   0.268    0.790    1.148    1.324    1.509
## beta_4                        1.694   0.335    1.019    1.474    1.707    1.926
## beta_5                        0.716   0.244    0.241    0.556    0.712    0.881
## beta_6                       -0.002   0.002   -0.007   -0.004   -0.002    0.000
## beta_7                        0.001   0.001   -0.002    0.000    0.001    0.002
## beta_8                        0.007   0.002    0.004    0.006    0.007    0.008
## beta_9                        0.000   0.001   -0.002   -0.001    0.000    0.001
## r                             4.240   0.265    3.743    4.045    4.241    4.419
## survival_non_censored[1]     84.868   0.000   84.868   84.868   84.868   84.868
## survival_non_censored[2]     65.947   0.000   65.947   65.947   65.947   65.947
## survival_non_censored[3]     80.857   0.000   80.857   80.857   80.857   80.857
## survival_non_censored[4]     81.517   0.000   81.517   81.517   81.517   81.517
## survival_non_censored[5]     82.601   0.000   82.601   82.601   82.601   82.601
## survival_non_censored[6]     81.695   0.000   81.695   81.695   81.695   81.695
## survival_non_censored[7]     67.168   0.000   67.168   67.168   67.168   67.168
## survival_non_censored[8]     79.214   0.000   79.214   79.214   79.214   79.214
## survival_non_censored[9]     93.380   0.000   93.380   93.380   93.380   93.380
## survival_non_censored[10]    85.736   0.000   85.736   85.736   85.736   85.736
## survival_non_censored[11]    80.698   0.000   80.698   80.698   80.698   80.698
## survival_non_censored[12]    69.027   0.000   69.027   69.027   69.027   69.027
## survival_non_censored[13]    68.468   0.000   68.468   68.468   68.468   68.468
## survival_non_censored[14]    81.013   0.000   81.013   81.013   81.013   81.013
## survival_non_censored[15]    81.676   0.000   81.676   81.676   81.676   81.676
## survival_non_censored[16]    68.893   0.000   68.893   68.893   68.893   68.893
## survival_non_censored[17]    75.907   0.000   75.907   75.907   75.907   75.907
## survival_non_censored[18]    83.088   0.000   83.088   83.088   83.088   83.088
## survival_non_censored[19]    87.830   0.000   87.830   87.830   87.830   87.830
## survival_non_censored[20]    81.049   0.000   81.049   81.049   81.049   81.049
## survival_non_censored[21]    68.816   0.000   68.816   68.816   68.816   68.816
## survival_non_censored[22]    71.652   0.000   71.652   71.652   71.652   71.652
## survival_non_censored[23]    85.541   0.000   85.541   85.541   85.541   85.541
## survival_non_censored[24]    80.780   0.000   80.780   80.780   80.780   80.780
## survival_non_censored[25]    78.242   0.000   78.242   78.242   78.242   78.242
## survival_non_censored[26]    86.026   0.000   86.026   86.026   86.026   86.026
## survival_non_censored[27]    69.864   0.000   69.864   69.864   69.864   69.864
## survival_non_censored[28]    68.268   0.000   68.268   68.268   68.268   68.268
## survival_non_censored[29]    80.674   0.000   80.674   80.674   80.674   80.674
## survival_non_censored[30]    76.315   0.000   76.315   76.315   76.315   76.315
## survival_non_censored[31]    69.478   0.000   69.478   69.478   69.478   69.478
## survival_non_censored[32]    70.366   0.000   70.366   70.366   70.366   70.366
## survival_non_censored[33]    69.862   0.000   69.862   69.862   69.862   69.862
## survival_non_censored[34]    69.021   0.000   69.021   69.021   69.021   69.021
## survival_non_censored[35]    72.446   0.000   72.446   72.446   72.446   72.446
## survival_non_censored[36]    56.676   0.000   56.676   56.676   56.676   56.676
## survival_non_censored[37]    69.147   0.000   69.147   69.147   69.147   69.147
## survival_non_censored[38]    68.704   0.000   68.704   68.704   68.704   68.704
## survival_non_censored[39]    83.255   0.000   83.255   83.255   83.255   83.255
## survival_non_censored[40]    68.287   0.000   68.287   68.287   68.287   68.287
## survival_non_censored[41]    66.691   0.000   66.691   66.691   66.691   66.691
## survival_non_censored[42]    83.135   0.000   83.135   83.135   83.135   83.135
## survival_non_censored[43]    53.985   0.000   53.985   53.985   53.985   53.985
## survival_non_censored[44]    67.529   0.000   67.529   67.529   67.529   67.529
## survival_non_censored[45]    81.695   0.000   81.695   81.695   81.695   81.695
## survival_non_censored[46]    56.331   0.000   56.331   56.331   56.331   56.331
## survival_non_censored[47]    64.534   0.000   64.534   64.534   64.534   64.534
## survival_non_censored[48]    45.971   0.000   45.971   45.971   45.971   45.971
## survival_non_censored[49]    69.213   0.000   69.213   69.213   69.213   69.213
## survival_non_censored[50]    64.386   0.000   64.386   64.386   64.386   64.386
## survival_non_censored[51]    72.624   0.000   72.624   72.624   72.624   72.624
## survival_non_censored[52]    60.564   0.000   60.564   60.564   60.564   60.564
## survival_non_censored[53]    70.062   0.000   70.062   70.062   70.062   70.062
## survival_non_censored[54]    54.418   0.000   54.418   54.418   54.418   54.418
## survival_non_censored[55]    58.823   0.000   58.823   58.823   58.823   58.823
## survival_non_censored[56]    79.595   0.000   79.595   79.595   79.595   79.595
## survival_non_censored[57]    57.355   0.000   57.355   57.355   57.355   57.355
## survival_non_censored[58]    64.142   0.000   64.142   64.142   64.142   64.142
## survival_non_censored[59]    62.133   0.000   62.133   62.133   62.133   62.133
## survival_non_censored[60]    91.135   0.000   91.135   91.135   91.135   91.135
## survival_non_censored[61]    67.844   0.000   67.844   67.844   67.844   67.844
## survival_non_censored[62]    67.808   0.000   67.808   67.808   67.808   67.808
## survival_non_censored[63]    90.675   0.000   90.675   90.675   90.675   90.675
## survival_non_censored[64]    83.222   0.000   83.222   83.222   83.222   83.222
## survival_non_censored[65]    85.284   0.000   85.284   85.284   85.284   85.284
## survival_non_censored[66]    73.180   0.000   73.180   73.180   73.180   73.180
## survival_non_censored[67]    80.619   0.000   80.619   80.619   80.619   80.619
## survival_non_censored[68]    78.231   0.000   78.231   78.231   78.231   78.231
## survival_non_censored[69]    79.630   0.000   79.630   79.630   79.630   79.630
## survival_non_censored[70]    68.145   0.000   68.145   68.145   68.145   68.145
## survival_non_censored[71]    71.806   0.000   71.806   71.806   71.806   71.806
## survival_non_censored[72]    53.615   0.000   53.615   53.615   53.615   53.615
## survival_non_censored[73]    65.618   0.000   65.618   65.618   65.618   65.618
## survival_non_censored[74]    74.163   0.000   74.163   74.163   74.163   74.163
## survival_non_censored[75]    64.873   0.000   64.873   64.873   64.873   64.873
## survival_non_censored[76]    77.106   0.000   77.106   77.106   77.106   77.106
## survival_non_censored[77]    56.170   0.000   56.170   56.170   56.170   56.170
## survival_non_censored[78]    66.585   0.000   66.585   66.585   66.585   66.585
## survival_non_censored[79]    63.239   0.000   63.239   63.239   63.239   63.239
## survival_non_censored[80]    70.289   0.000   70.289   70.289   70.289   70.289
## survival_non_censored[81]    49.834   0.000   49.834   49.834   49.834   49.834
## survival_non_censored[82]    57.120   0.000   57.120   57.120   57.120   57.120
## survival_non_censored[83]    71.266   0.000   71.266   71.266   71.266   71.266
## survival_non_censored[84]    67.559   0.000   67.559   67.559   67.559   67.559
## survival_non_censored[85]    58.623   0.000   58.623   58.623   58.623   58.623
## survival_non_censored[86]    60.192   0.000   60.192   60.192   60.192   60.192
## survival_non_censored[87]    72.474   0.000   72.474   72.474   72.474   72.474
## survival_non_censored[88]    67.097   0.000   67.097   67.097   67.097   67.097
## survival_non_censored[89]    57.744   0.000   57.744   57.744   57.744   57.744
## survival_non_censored[90]    60.504   0.000   60.504   60.504   60.504   60.504
## survival_non_censored[91]    90.193   0.000   90.193   90.193   90.193   90.193
## survival_non_censored[92]    63.195   0.000   63.195   63.195   63.195   63.195
## survival_non_censored[93]    88.632   0.000   88.632   88.632   88.632   88.632
## survival_non_censored[94]    78.450   0.000   78.450   78.450   78.450   78.450
## survival_non_censored[95]    46.483   0.000   46.483   46.483   46.483   46.483
## survival_non_censored[96]    64.405   0.000   64.405   64.405   64.405   64.405
## survival_non_censored[97]    81.227   0.000   81.227   81.227   81.227   81.227
## survival_non_censored[98]    93.451   0.000   93.451   93.451   93.451   93.451
## survival_non_censored[99]    93.328   0.000   93.328   93.328   93.328   93.328
## survival_non_censored[100]   94.467   0.000   94.467   94.467   94.467   94.467
## survival_non_censored[101]   69.673   0.000   69.673   69.673   69.673   69.673
## survival_non_censored[102]   64.277   0.000   64.277   64.277   64.277   64.277
## survival_non_censored[103]   46.782   0.000   46.782   46.782   46.782   46.782
## survival_non_censored[104]   35.877   0.000   35.877   35.877   35.877   35.877
## survival_non_censored[105]   36.235   0.000   36.235   36.235   36.235   36.235
## survival_non_censored[106]   28.476   0.000   28.476   28.476   28.476   28.476
## survival_non_censored[107]   39.751   0.000   39.751   39.751   39.751   39.751
## survival_non_censored[108]   34.856   0.000   34.856   34.856   34.856   34.856
## survival_non_censored[109]   29.478   0.000   29.478   29.478   29.478   29.478
## survival_non_censored[110]   60.186   0.000   60.186   60.186   60.186   60.186
## survival_non_censored[111]   35.337   0.000   35.337   35.337   35.337   35.337
## survival_non_censored[112]   56.956   0.000   56.956   56.956   56.956   56.956
## survival_non_censored[113]   38.081   0.000   38.081   38.081   38.081   38.081
## survival_non_censored[114]   21.769   0.000   21.769   21.769   21.769   21.769
## survival_non_censored[115]   33.243   0.000   33.243   33.243   33.243   33.243
## survival_non_censored[116]   22.897   0.000   22.897   22.897   22.897   22.897
## survival_non_censored[117]   68.624   0.000   68.624   68.624   68.624   68.624
## survival_non_censored[118]   56.816   0.000   56.816   56.816   56.816   56.816
## survival_non_censored[119]   87.203   0.000   87.203   87.203   87.203   87.203
## survival_non_censored[120]   59.800   0.000   59.800   59.800   59.800   59.800
## survival_non_censored[121]   67.444   0.000   67.444   67.444   67.444   67.444
## survival_non_censored[122]   30.100   0.000   30.100   30.100   30.100   30.100
## survival_non_censored[123]   18.710   0.000   18.710   18.710   18.710   18.710
## survival_non_censored[124]   37.251   0.000   37.251   37.251   37.251   37.251
## survival_non_censored[125]   61.689   0.000   61.689   61.689   61.689   61.689
## survival_non_censored[126]   82.998   0.000   82.998   82.998   82.998   82.998
## survival_non_censored[127]   66.998   0.000   66.998   66.998   66.998   66.998
## survival_non_censored[128]   44.999   0.000   44.999   44.999   44.999   44.999
## survival_non_censored[129]   27.907   0.000   27.907   27.907   27.907   27.907
## survival_non_censored[130]   64.999   0.000   64.999   64.999   64.999   64.999
## survival_non_censored[131]   23.707   0.000   23.707   23.707   23.707   23.707
## survival_non_censored[132]   49.002   0.000   49.002   49.002   49.002   49.002
## survival_non_censored[133]   45.996   0.000   45.996   45.996   45.996   45.996
## survival_non_censored[134]    9.259   0.000    9.259    9.259    9.259    9.259
## survival_non_censored[135]   21.506   0.000   21.506   21.506   21.506   21.506
## survival_non_censored[136]   17.248   0.000   17.248   17.248   17.248   17.248
## survival_non_censored[137]   18.242   0.000   18.242   18.242   18.242   18.242
## survival_non_censored[138]   57.843   0.000   57.843   57.843   57.843   57.843
## survival_non_censored[139]   40.241   0.000   40.241   40.241   40.241   40.241
## survival_non_censored[140]   51.652   0.000   51.652   51.652   51.652   51.652
## survival_non_censored[141]   77.352   0.000   77.352   77.352   77.352   77.352
## survival_non_censored[142]   56.331   0.000   56.331   56.331   56.331   56.331
## survival_non_censored[143]   27.302   0.000   27.302   27.302   27.302   27.302
## survival_non_censored[144]   51.526   0.000   51.526   51.526   51.526   51.526
## survival_non_censored[145]   58.300   0.000   58.300   58.300   58.300   58.300
## survival_non_censored[146]   63.493   0.000   63.493   63.493   63.493   63.493
## survival_non_censored[147]   62.667   0.000   62.667   62.667   62.667   62.667
## survival_non_censored[148]   75.639   0.000   75.639   75.639   75.639   75.639
## survival_non_censored[149]   45.818   0.000   45.818   45.818   45.818   45.818
## survival_non_censored[150]   84.203   0.000   84.203   84.203   84.203   84.203
## survival_non_censored[151]   72.903   0.000   72.903   72.903   72.903   72.903
## survival_non_censored[152]   21.528   0.000   21.528   21.528   21.528   21.528
## survival_non_censored[153]   47.220   0.000   47.220   47.220   47.220   47.220
## survival_non_censored[154]   78.209   0.000   78.209   78.209   78.209   78.209
## survival_non_censored[155]   34.237   0.000   34.237   34.237   34.237   34.237
## survival_non_censored[156]   35.318   0.000   35.318   35.318   35.318   35.318
## survival_non_censored[157]   30.300   0.000   30.300   30.300   30.300   30.300
## survival_non_censored[158]   21.380   0.000   21.380   21.380   21.380   21.380
## survival_non_censored[159]   73.248   0.000   73.248   73.248   73.248   73.248
## survival_non_censored[160]   21.363   0.000   21.363   21.363   21.363   21.363
## survival_non_censored[161]   69.216   0.000   69.216   69.216   69.216   69.216
## survival_non_censored[162]   45.936   0.000   45.936   45.936   45.936   45.936
## survival_non_censored[163]   35.526   0.000   35.526   35.526   35.526   35.526
## survival_non_censored[164]   59.734   0.000   59.734   59.734   59.734   59.734
## survival_non_censored[165]   47.316   0.000   47.316   47.316   47.316   47.316
## survival_non_censored[166]   87.693   0.000   87.693   87.693   87.693   87.693
## deviance                   1423.309   4.725 1416.048 1419.895 1422.654 1426.055
##                               97.5%  Rhat n.eff
## age_Dalai_predictive        118.774 1.002  1100
## age_Francis_predictive      119.424 1.001  3000
## age_Naruhito_predictive     119.328 1.001  2900
## age_Obama_predictive        119.270 1.001  3000
## beta_0                      -19.587 1.051    53
## beta_1                        0.000 1.002  1100
## beta_2                        1.416 1.013   160
## beta_3                        1.848 1.009   250
## beta_4                        2.331 1.022   110
## beta_5                        1.201 1.012   170
## beta_6                        0.003 1.003   740
## beta_7                        0.003 1.001  2400
## beta_8                        0.010 1.003   890
## beta_9                        0.003 1.002  1700
## r                             4.762 1.027   110
## survival_non_censored[1]     84.868 1.000     1
## survival_non_censored[2]     65.947 1.000     1
## survival_non_censored[3]     80.857 1.000     1
## survival_non_censored[4]     81.517 1.000     1
## survival_non_censored[5]     82.601 1.000     1
## survival_non_censored[6]     81.695 1.000     1
## survival_non_censored[7]     67.168 1.000     1
## survival_non_censored[8]     79.214 1.000     1
## survival_non_censored[9]     93.380 1.000     1
## survival_non_censored[10]    85.736 1.000     1
## survival_non_censored[11]    80.698 1.000     1
## survival_non_censored[12]    69.027 1.000     1
## survival_non_censored[13]    68.468 1.000     1
## survival_non_censored[14]    81.013 1.000     1
## survival_non_censored[15]    81.676 1.000     1
## survival_non_censored[16]    68.893 1.000     1
## survival_non_censored[17]    75.907 1.000     1
## survival_non_censored[18]    83.088 1.000     1
## survival_non_censored[19]    87.830 1.000     1
## survival_non_censored[20]    81.049 1.000     1
## survival_non_censored[21]    68.816 1.000     1
## survival_non_censored[22]    71.652 1.000     1
## survival_non_censored[23]    85.541 1.000     1
## survival_non_censored[24]    80.780 1.000     1
## survival_non_censored[25]    78.242 1.000     1
## survival_non_censored[26]    86.026 1.000     1
## survival_non_censored[27]    69.864 1.000     1
## survival_non_censored[28]    68.268 1.000     1
## survival_non_censored[29]    80.674 1.000     1
## survival_non_censored[30]    76.315 1.000     1
## survival_non_censored[31]    69.478 1.000     1
## survival_non_censored[32]    70.366 1.000     1
## survival_non_censored[33]    69.862 1.000     1
## survival_non_censored[34]    69.021 1.000     1
## survival_non_censored[35]    72.446 1.000     1
## survival_non_censored[36]    56.676 1.000     1
## survival_non_censored[37]    69.147 1.000     1
## survival_non_censored[38]    68.704 1.000     1
## survival_non_censored[39]    83.255 1.000     1
## survival_non_censored[40]    68.287 1.000     1
## survival_non_censored[41]    66.691 1.000     1
## survival_non_censored[42]    83.135 1.000     1
## survival_non_censored[43]    53.985 1.000     1
## survival_non_censored[44]    67.529 1.000     1
## survival_non_censored[45]    81.695 1.000     1
## survival_non_censored[46]    56.331 1.000     1
## survival_non_censored[47]    64.534 1.000     1
## survival_non_censored[48]    45.971 1.000     1
## survival_non_censored[49]    69.213 1.000     1
## survival_non_censored[50]    64.386 1.000     1
## survival_non_censored[51]    72.624 1.000     1
## survival_non_censored[52]    60.564 1.000     1
## survival_non_censored[53]    70.062 1.000     1
## survival_non_censored[54]    54.418 1.000     1
## survival_non_censored[55]    58.823 1.000     1
## survival_non_censored[56]    79.595 1.000     1
## survival_non_censored[57]    57.355 1.000     1
## survival_non_censored[58]    64.142 1.000     1
## survival_non_censored[59]    62.133 1.000     1
## survival_non_censored[60]    91.135 1.000     1
## survival_non_censored[61]    67.844 1.000     1
## survival_non_censored[62]    67.808 1.000     1
## survival_non_censored[63]    90.675 1.000     1
## survival_non_censored[64]    83.222 1.000     1
## survival_non_censored[65]    85.284 1.000     1
## survival_non_censored[66]    73.180 1.000     1
## survival_non_censored[67]    80.619 1.000     1
## survival_non_censored[68]    78.231 1.000     1
## survival_non_censored[69]    79.630 1.000     1
## survival_non_censored[70]    68.145 1.000     1
## survival_non_censored[71]    71.806 1.000     1
## survival_non_censored[72]    53.615 1.000     1
## survival_non_censored[73]    65.618 1.000     1
## survival_non_censored[74]    74.163 1.000     1
## survival_non_censored[75]    64.873 1.000     1
## survival_non_censored[76]    77.106 1.000     1
## survival_non_censored[77]    56.170 1.000     1
## survival_non_censored[78]    66.585 1.000     1
## survival_non_censored[79]    63.239 1.000     1
## survival_non_censored[80]    70.289 1.000     1
## survival_non_censored[81]    49.834 1.000     1
## survival_non_censored[82]    57.120 1.000     1
## survival_non_censored[83]    71.266 1.000     1
## survival_non_censored[84]    67.559 1.000     1
## survival_non_censored[85]    58.623 1.000     1
## survival_non_censored[86]    60.192 1.000     1
## survival_non_censored[87]    72.474 1.000     1
## survival_non_censored[88]    67.097 1.000     1
## survival_non_censored[89]    57.744 1.000     1
## survival_non_censored[90]    60.504 1.000     1
## survival_non_censored[91]    90.193 1.000     1
## survival_non_censored[92]    63.195 1.000     1
## survival_non_censored[93]    88.632 1.000     1
## survival_non_censored[94]    78.450 1.000     1
## survival_non_censored[95]    46.483 1.000     1
## survival_non_censored[96]    64.405 1.000     1
## survival_non_censored[97]    81.227 1.000     1
## survival_non_censored[98]    93.451 1.000     1
## survival_non_censored[99]    93.328 1.000     1
## survival_non_censored[100]   94.467 1.000     1
## survival_non_censored[101]   69.673 1.000     1
## survival_non_censored[102]   64.277 1.000     1
## survival_non_censored[103]   46.782 1.000     1
## survival_non_censored[104]   35.877 1.000     1
## survival_non_censored[105]   36.235 1.000     1
## survival_non_censored[106]   28.476 1.000     1
## survival_non_censored[107]   39.751 1.000     1
## survival_non_censored[108]   34.856 1.000     1
## survival_non_censored[109]   29.478 1.000     1
## survival_non_censored[110]   60.186 1.000     1
## survival_non_censored[111]   35.337 1.000     1
## survival_non_censored[112]   56.956 1.000     1
## survival_non_censored[113]   38.081 1.000     1
## survival_non_censored[114]   21.769 1.000     1
## survival_non_censored[115]   33.243 1.000     1
## survival_non_censored[116]   22.897 1.000     1
## survival_non_censored[117]   68.624 1.000     1
## survival_non_censored[118]   56.816 1.000     1
## survival_non_censored[119]   87.203 1.000     1
## survival_non_censored[120]   59.800 1.000     1
## survival_non_censored[121]   67.444 1.000     1
## survival_non_censored[122]   30.100 1.000     1
## survival_non_censored[123]   18.710 1.000     1
## survival_non_censored[124]   37.251 1.000     1
## survival_non_censored[125]   61.689 1.000     1
## survival_non_censored[126]   82.998 1.000     1
## survival_non_censored[127]   66.998 1.000     1
## survival_non_censored[128]   44.999 1.000     1
## survival_non_censored[129]   27.907 1.000     1
## survival_non_censored[130]   64.999 1.000     1
## survival_non_censored[131]   23.707 1.000     1
## survival_non_censored[132]   49.002 1.000     1
## survival_non_censored[133]   45.996 1.000     1
## survival_non_censored[134]    9.259 1.000     1
## survival_non_censored[135]   21.506 1.000     1
## survival_non_censored[136]   17.248 1.000     1
## survival_non_censored[137]   18.242 1.000     1
## survival_non_censored[138]   57.843 1.000     1
## survival_non_censored[139]   40.241 1.000     1
## survival_non_censored[140]   51.652 1.000     1
## survival_non_censored[141]   77.352 1.000     1
## survival_non_censored[142]   56.331 1.000     1
## survival_non_censored[143]   27.302 1.000     1
## survival_non_censored[144]   51.526 1.000     1
## survival_non_censored[145]   58.300 1.000     1
## survival_non_censored[146]   63.493 1.000     1
## survival_non_censored[147]   62.667 1.000     1
## survival_non_censored[148]   75.639 1.000     1
## survival_non_censored[149]   45.818 1.000     1
## survival_non_censored[150]   84.203 1.000     1
## survival_non_censored[151]   72.903 1.000     1
## survival_non_censored[152]   21.528 1.000     1
## survival_non_censored[153]   47.220 1.000     1
## survival_non_censored[154]   78.209 1.000     1
## survival_non_censored[155]   34.237 1.000     1
## survival_non_censored[156]   35.318 1.000     1
## survival_non_censored[157]   30.300 1.000     1
## survival_non_censored[158]   21.380 1.000     1
## survival_non_censored[159]   73.248 1.000     1
## survival_non_censored[160]   21.363 1.000     1
## survival_non_censored[161]   69.216 1.000     1
## survival_non_censored[162]   45.936 1.000     1
## survival_non_censored[163]   35.526 1.000     1
## survival_non_censored[164]   59.734 1.000     1
## survival_non_censored[165]   47.316 1.000     1
## survival_non_censored[166]   87.693 1.000     1
## deviance                   1433.614 1.010   250
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## 
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 11.1 and DIC = 1434.4
## DIC is an estimate of expected predictive error (lower deviance is better).

library(coda)
parameters = c("beta_0", "beta_1", "beta_2", "beta_3", "beta_4", "beta_5", 
               "beta_6", "beta_7",  "beta_8", "beta_9", "r")

res = data.frame(Bayesian_Popes_second_approach$BUGSoutput$sims.matrix)[,parameters]
colnames(res) <- gsub("^.*\\.","", colnames(res) )

# lapply(seq(1,length(res)), function(i) {
#   traceplot(as.mcmc(res[,i]), smooth = FALSE, type = "l",
#             xlab = "Iterations", ylab = colnames(res)[i],
#             main = paste("Traceplot of ", colnames(res)[i]))
# })

tps <- function(var){
  ggplot(res, aes_(y=as.name(var), x=seq(1,nrow(res)))) +
    geom_line() +
    labs(title=paste("Traceplot of ", as.name(var)),
         x ="Iterations", y = as.name(var))
}
lapply(names(res), tps)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

library(ggplot2)
res_lag1 = lapply(seq(1:length(res)), function(i) {
  lres = lag(res[,i],1)
  plot(y=res[,i], x= lres, 
       xlab = paste0(colnames(res)[i], "-1"),
       ylab = paste0(colnames(res)[i]),
       main = paste("Lag-1 Scatter Plot of", colnames(res)[i]))
})

lapply(seq(1,length(res)), function(i) { 
  acf(res[,i], xlab = "Lag", ylab = "ACF", 
            main = paste("ACF Plot of ", colnames(res)[i]))
})

## [[1]]
## 
## Autocorrelations of series 'res[, i]', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.017  0.037 -0.022  0.008 -0.004  0.022  0.022 -0.018  0.001  0.039 
##     11     12     13     14     15     16     17     18     19     20     21 
## -0.003 -0.004  0.007 -0.018 -0.012  0.003  0.004 -0.001  0.003 -0.027 -0.009 
##     22     23     24     25     26     27     28     29     30     31     32 
##  0.016 -0.020 -0.004 -0.012 -0.020  0.028  0.010 -0.013 -0.016  0.008  0.010 
##     33     34 
##  0.000  0.021 
## 
## [[2]]
## 
## Autocorrelations of series 'res[, i]', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.001  0.039 -0.027  0.002  0.016 -0.002  0.006  0.029 -0.014 -0.002 
##     11     12     13     14     15     16     17     18     19     20     21 
## -0.006  0.020  0.015 -0.005  0.011 -0.021  0.003 -0.016  0.024 -0.003 -0.024 
##     22     23     24     25     26     27     28     29     30     31     32 
## -0.015  0.005  0.051 -0.022 -0.007  0.027  0.014  0.007 -0.039 -0.019 -0.034 
##     33     34 
## -0.016  0.006 
## 
## [[3]]
## 
## Autocorrelations of series 'res[, i]', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000 -0.005  0.008  0.020  0.003  0.017  0.039  0.000 -0.028  0.011  0.008 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.008  0.002 -0.019 -0.002  0.005  0.016  0.016 -0.019  0.021 -0.026  0.008 
##     22     23     24     25     26     27     28     29     30     31     32 
## -0.001  0.004  0.037  0.007  0.009  0.015  0.004  0.014  0.001  0.015  0.006 
##     33     34 
##  0.020  0.002 
## 
## [[4]]
## 
## Autocorrelations of series 'res[, i]', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000 -0.036  0.029 -0.008  0.016  0.035  0.002  0.009  0.027 -0.019  0.028 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.007  0.005 -0.005  0.005  0.002  0.003 -0.023  0.010 -0.004  0.020  0.016 
##     22     23     24     25     26     27     28     29     30     31     32 
##  0.015 -0.015  0.009 -0.003 -0.034  0.000  0.021 -0.028 -0.025 -0.023  0.017 
##     33     34 
##  0.007  0.008 
## 
## [[5]]
## 
## Autocorrelations of series 'res[, i]', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.024 -0.001  0.010  0.005 -0.030  0.005  0.016  0.007 -0.010  0.003 
##     11     12     13     14     15     16     17     18     19     20     21 
## -0.006 -0.032 -0.009  0.009 -0.008 -0.020 -0.016 -0.027 -0.017 -0.006 -0.004 
##     22     23     24     25     26     27     28     29     30     31     32 
## -0.011  0.004  0.039  0.031 -0.022 -0.018 -0.002  0.011  0.002  0.005 -0.004 
##     33     34 
## -0.029  0.026 
## 
## [[6]]
## 
## Autocorrelations of series 'res[, i]', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.002 -0.020 -0.026  0.005 -0.005  0.009  0.014  0.016  0.007  0.016 
##     11     12     13     14     15     16     17     18     19     20     21 
## -0.017 -0.028  0.016  0.002  0.001 -0.030 -0.002 -0.019  0.000 -0.006  0.012 
##     22     23     24     25     26     27     28     29     30     31     32 
##  0.009 -0.038  0.004  0.008 -0.019  0.007 -0.035  0.024 -0.012 -0.010 -0.002 
##     33     34 
##  0.025  0.021 
## 
## [[7]]
## 
## Autocorrelations of series 'res[, i]', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000 -0.026  0.012  0.008  0.023 -0.010  0.014 -0.009 -0.030 -0.011  0.007 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.015  0.001 -0.021  0.028  0.007 -0.010  0.013 -0.013  0.034 -0.034  0.011 
##     22     23     24     25     26     27     28     29     30     31     32 
##  0.012  0.016  0.026  0.010 -0.010  0.000  0.031  0.024 -0.004  0.015  0.002 
##     33     34 
##  0.016  0.038 
## 
## [[8]]
## 
## Autocorrelations of series 'res[, i]', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000 -0.015  0.013 -0.003 -0.008 -0.014 -0.037 -0.013  0.005  0.021  0.001 
##     11     12     13     14     15     16     17     18     19     20     21 
## -0.008  0.043  0.009  0.019 -0.017  0.023 -0.021 -0.002 -0.013 -0.001 -0.035 
##     22     23     24     25     26     27     28     29     30     31     32 
## -0.027  0.012  0.034  0.020  0.010 -0.006  0.012  0.002 -0.008  0.009 -0.022 
##     33     34 
## -0.013  0.011 
## 
## [[9]]
## 
## Autocorrelations of series 'res[, i]', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000 -0.023  0.039 -0.012 -0.005 -0.022 -0.038 -0.003  0.022 -0.026 -0.015 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.001  0.008  0.012  0.015  0.003  0.015 -0.031 -0.003  0.018  0.006  0.021 
##     22     23     24     25     26     27     28     29     30     31     32 
##  0.000 -0.036  0.018  0.008 -0.025  0.031  0.018  0.006 -0.005 -0.033  0.007 
##     33     34 
##  0.007 -0.015 
## 
## [[10]]
## 
## Autocorrelations of series 'res[, i]', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.022  0.022  0.006 -0.001  0.019  0.002  0.035  0.016 -0.012 -0.006 
##     11     12     13     14     15     16     17     18     19     20     21 
## -0.016  0.021 -0.006 -0.016  0.009 -0.003 -0.011  0.012  0.023  0.006 -0.053 
##     22     23     24     25     26     27     28     29     30     31     32 
##  0.028 -0.002  0.022 -0.029 -0.017  0.011 -0.016  0.012 -0.008 -0.043 -0.023 
##     33     34 
## -0.026  0.014 
## 
## [[11]]
## 
## Autocorrelations of series 'res[, i]', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.014  0.011 -0.013  0.008  0.015  0.003  0.013 -0.014  0.002  0.027 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.007 -0.013  0.014 -0.009 -0.010  0.009  0.012  0.028  0.005  0.009 -0.020 
##     22     23     24     25     26     27     28     29     30     31     32 
## -0.011 -0.012 -0.019 -0.018 -0.020  0.024 -0.006 -0.003  0.006  0.012  0.003 
##     33     34 
## -0.016  0.002

people = c("age_Francis_predictive", "age_Obama_predictive", "age_Dalai_predictive",
"age_Naruhito_predictive")
four_pred = data.frame(Bayesian_Popes_second_approach$BUGSoutput$sims.matrix)[,people]

print(paste("P(Dalai Lifespan > Francis Lifespan)=", mean(four_pred$age_Dalai_predictive>four_pred$age_Francis_predictive)))

## [1] "P(Dalai Lifespan > Francis Lifespan)= 0.383666666666667"

print(paste("P(Obama Lifespan > Naruhito Lifespan)=", mean(four_pred$age_Obama_predictive>four_pred$age_Naruhito_predictive)))

## [1] "P(Obama Lifespan > Naruhito Lifespan)= 0.501"

Posterior Predictive Plots (4 separate)

ppd_plot <- function(i) {
  name = gsub(".*[_]([^.]+)[_].*", "\\1", colnames(four_pred)[i])
  minval = round(min(four_pred[,i]), 2)
  medval = round(median(four_pred[,i]), 2)
  val95 = round(quantile(four_pred[,i], 0.95),2)
  maxval = round(max(four_pred[,i]), 2)
ggplot(four_pred, aes(x=four_pred[,i])) + geom_histogram(aes(y=..density..), 
                                                   binwidth = 0.8,
                                                   boundary = minval) +
  geom_vline(xintercept = minval, colour = "black")+
  geom_vline(xintercept = medval, colour = "black") +
  geom_vline(xintercept = val95, colour = "black")+
  geom_hline(yintercept = 0, colour = "black") +
  xlab("Lifespan (years)") +
  ylab("Predictive Probability density") +
  ggtitle(paste("Predictive Posterior Distribution for Lifespan of", name)) +
  ylim(c(0, 0.07))+
  annotate("text", x = minval+0.15*(medval-minval), y = 0.06, label = paste("Current age", minval, "years", sep = "\n"), size = 3)+
  annotate("text", x = (medval+minval)/2, y = 0.05, label = "50%", size = 3)+
  annotate("text", x = medval+0.15*(val95-medval), y = 0.06, label = paste("Median", medval,"years", sep = "\n"), size = 3)+
  annotate("text", x = (medval+val95)/2, y = 0.05, label = "45%", size = 3)+
  annotate("text", x = val95+2.2, y = 0.06, label = paste("95%", "quantile", val95, "years", sep = "\n"), size = 3)+
  annotate("text", x = val95+3, y = 0.05, label = "5%", size = 3)
}

lapply(seq(1,4), ppd_plot)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

Posterior Predictive Plots (2 pairs)

Francis and Dalai Lama

francis_dalai = data.frame(plifespans = c(four_pred$age_Francis_predictive, four_pred$age_Dalai_predictive), person = c(rep("Francis", nrow(four_pred)), rep("Dalai", nrow(four_pred))))

mins = francis_dalai %>%
  group_by(person) %>%
  summarise(min_vals = round(min(plifespans), 2))

## `summarise()` ungrouping output (override with `.groups` argument)

meds = francis_dalai %>%
  group_by(person) %>%
  summarise(med_vals = round(median(plifespans), 2))

## `summarise()` ungrouping output (override with `.groups` argument)

val95 = francis_dalai %>%
  group_by(person) %>%
  summarise(vals_95 = round(quantile(plifespans, 0.95), 2))

## `summarise()` ungrouping output (override with `.groups` argument)

francis_txt = paste(paste("Current age:", mins$min_vals[mins$person=="Francis"], "yrs"),
      paste("Median:", meds$med_vals[meds$person=="Francis"],"yrs"),
      paste("95% quantile:", val95$vals_95[val95$person=="Francis"], "yrs"),
      sep = "\n")

dalai_txt = paste(paste("Current age:", mins$min_vals[mins$person=="Dalai"], "yrs"),
      paste("Median:",  meds$med_vals[meds$person=="Dalai"],"yrs"),
      paste("95% quantile:", val95$vals_95[val95$person=="Dalai"], "yrs"),
      sep = "\n")

plot_txt <- data.frame(
  label = c(francis_txt, dalai_txt),
  person = c("Francis", "Dalai")
)

ggplot(francis_dalai, aes(x=plifespans)) + 
  geom_histogram(aes(y=..density..), binwidth = 0.8,boundary = round(min(francis_dalai$plifespans), 2)) +
  facet_wrap(~person) +
  geom_vline(data=mins, aes(xintercept=min_vals))+
  geom_vline(data=meds, aes(xintercept=med_vals))+
  geom_vline(data=val95, aes(xintercept=vals_95))+
  geom_hline(yintercept = 0, colour = "black") +
  xlab("Lifespan (yrs)") + ylab("Predictive Probability density") +
  ggtitle(paste("Posterior Predictive Lifespans (Pope Francis and 14th Dalai Lama)")) +
  ylim(c(0, 0.05)) + 
  geom_text(data = plot_txt, mapping = aes(x = 92.8, y = 0.048, label = label), size = 2.8)

Obama and Naruhito

obama_naru = data.frame(plifespans = c(four_pred$age_Obama_predictive, four_pred$age_Naruhito_predictive), person = c(rep("Obama", nrow(four_pred)), rep("Naruhito", nrow(four_pred))))

mins = obama_naru %>%
  group_by(person) %>%
  summarise(min_vals = round(min(plifespans), 2))

## `summarise()` ungrouping output (override with `.groups` argument)

meds = obama_naru %>%
  group_by(person) %>%
  summarise(med_vals = round(median(plifespans), 2))

## `summarise()` ungrouping output (override with `.groups` argument)

val95 = obama_naru %>%
  group_by(person) %>%
  summarise(vals_95 = round(quantile(plifespans, 0.95), 2))

## `summarise()` ungrouping output (override with `.groups` argument)

obama_txt = paste(paste("Current age:", mins$min_vals[mins$person=="Obama"], "yrs"),
      paste("Median:", meds$med_vals[meds$person=="Obama"],"yrs"),
      paste("95% quantile:", val95$vals_95[val95$person=="Obama"], "yrs"),
      sep = "\n")

naru_txt = paste(paste("Current age:", mins$min_vals[mins$person=="Naruhito"], "yrs"),
      paste("Median:",  meds$med_vals[meds$person=="Naruhito"],"yrs"),
      paste("95% quantile:", val95$vals_95[val95$person=="Naruhito"], "yrs"),
      sep = "\n")

plot_txt <- data.frame(
  label = c(obama_txt, naru_txt),
  person = c("Obama", "Naruhito")
)

ggplot(obama_naru, aes(x=plifespans)) + 
  geom_histogram(aes(y=..density..), binwidth = 0.8,boundary = round(min(obama_naru$plifespans), 2)) +
  facet_wrap(~person) +
  geom_vline(data=mins, aes(xintercept=min_vals))+
  geom_vline(data=meds, aes(xintercept=med_vals))+
  geom_vline(data=val95, aes(xintercept=vals_95))+
  geom_hline(yintercept = 0, colour = "black") +
  xlab("Lifespan (yrs)") + ylab("Predictive Probability density") +
  ggtitle(paste("Posterior Predictive Lifespans (President Obama and Emperor Naruhito)")) +
  ylim(c(0, 0.05)) + 
  geom_text(data = plot_txt, mapping = aes(x = 80, y = 0.04, label = label), size = 3)

vars = data.frame(Bayesian_Popes_second_approach$BUGSoutput$root.short)
output <- head(data.frame(Bayesian_Popes_second_approach$BUGSoutput$summary), 16)
output <- output %>%
  mutate("Variable" = head(vars,16),
         "Mean" = output$mean,
         "Standard Deviation" = output$sd,
         "2.5% Quantile" = output$X2.5.,
         "Median" = output$X50.,
         "97.5% Quantile" = output$X97.5.
         ) %>%
  select("Variable", "Mean", "Standard Deviation", "2.5% Quantile", "Median", "97.5% Quantile")
coef <- output %>%
  slice(5:16)
kable(coef)

Variable	Mean	Standard Deviation	2.5% Quantile	Median	97.5% Quantile
beta_0	-22.7763077	1.5962011	-25.9407203	-22.7453356	-19.5873392
beta_1	-0.0016375	0.0007812	-0.0031670	-0.0016413	-0.0001422
beta_2	0.6008876	0.4402710	-0.2822430	0.6097154	1.4157974
beta_3	1.3236079	0.2677447	0.7903557	1.3239796	1.8480258
beta_4	1.6942849	0.3354591	1.0192361	1.7074271	2.3307992
beta_5	0.7160662	0.2440568	0.2412665	0.7124606	1.2012485
beta_6	-0.0020120	0.0024374	-0.0067747	-0.0020312	0.0026784
beta_7	0.0006767	0.0013981	-0.0020156	0.0006623	0.0034650
beta_8	0.0068935	0.0017716	0.0036422	0.0068926	0.0103842
beta_9	0.0002039	0.0012350	-0.0021660	0.0001918	0.0025679
deviance	1423.3093231	4.7248506	1416.0482075	1422.6535327	1433.6141160
r	4.2402583	0.2647113	3.7428743	4.2410108	4.7623073

beta_0_est <- coef[1,2]
beta_1_est <- coef[2,2]
beta_2_est <- coef[3,2]
beta_3_est <- coef[4,2]
beta_4_est <- coef[5,2]
beta_5_est <- coef[6,2]
beta_6_est <- coef[7,2]
beta_7_est <- coef[8,2]
beta_8_est <- coef[9,2]
beta_9_est <- coef[10,2]

no_living <- leaders_data %>%
  filter(Censored == 0)

estimates <- no_living %>%
  mutate(Birth.Year.Cent = Birth.Year - mean(leaders_nopred$Birth.Year)) %>%
  mutate(ChinaEmp =case_when(Type == "ChinaEmp" ~ 1, TRUE ~ 0),
         DalaiLama =case_when(Type == "DalaiLama" ~ 1, TRUE ~ 0),
         JapanEmp =case_when(Type == "JapanEmp" ~ 1, TRUE ~ 0),
         Pope =case_when(Type == "Pope" ~ 1, TRUE ~ 0),
         UsPres =case_when(Type == "UsPres" ~ 1, TRUE ~ 0)) %>%
  mutate(BYUsPres = Birth.Year.Cent*UsPres) %>%
  mutate(BYChinaEmp = Birth.Year.Cent*ChinaEmp) %>%
  mutate(BYDalaiLama = Birth.Year.Cent*DalaiLama) %>%
  mutate(BYJapanEmp = Birth.Year.Cent*JapanEmp) %>%
  mutate(intercept = 1) %>%
  mutate(predicted = exp(beta_0_est + beta_1_est*(Birth.Year.Cent) + beta_2_est*(UsPres) + beta_3_est*(ChinaEmp) + beta_4_est*(DalaiLama) + beta_5_est*(JapanEmp) + beta_6_est*(UsPres)*(Birth.Year.Cent) + beta_7_est*(ChinaEmp)*(Birth.Year.Cent) + beta_8_est*(DalaiLama)*(Birth.Year.Cent) + beta_9_est*(JapanEmp)*(Birth.Year.Cent)))


# check using matrices -- Ts in the 200s
data.mat = as.matrix(estimates[,c("intercept", "Birth.Year.Cent", "UsPres", "ChinaEmp", "DalaiLama", "JapanEmp","BYUsPres", "BYChinaEmp", "BYDalaiLama", "BYJapanEmp")], nrow = nrow(estimates), byrow = TRUE)

r = Bayesian_Popes_second_approach$BUGSoutput$summary[,1][16]
betas = Bayesian_Popes_second_approach$BUGSoutput$summary[,1][5:14]
beta.mat = as.matrix(betas, nrow = p, byrow = TRUE)
logmus = data.mat %*% beta.mat

min.age = min(leaders_nopred$Age.Event)

install.packages("TruncatedDistributions", repos="http://R-Forge.R-project.org")

## Installing package into '/home/rstudio-user/R/x86_64-pc-linux-gnu-library/4.0'
## (as 'lib' is unspecified)

library(TruncatedDistributions)
Ts = rtweibull(166, r, (1/exp(logmus))^(1/r), min.age, 120) 
# dweibull takes a diff parameterization of Weibull
# please check the parameterization see ?dweibull

Posterior Predictive Checks

non_censored_prediction <- filter(leaders_data, Censored == 0)
K <- 10
ysims <- matrix(nrow = nrow(res), ncol = K)
for(k in 1:K){
  for(i in 1:nrow(non_censored_prediction)){
    w <- sample(nrow(res), 1)
    samp <- res[w,]
    val <- non_censored_prediction[i, ]
    beta_temp <- samp$beta_0 + samp$beta_1*val$Birth.Year + samp$beta_2*as.integer(val$TypeUsPres) + samp$beta_3*as.integer(val$TypeChinaEmp) + samp$beta_4*as.integer(val$TypeDalaiLama) + samp$beta_5*as.integer(val$TypeJapanEmp) + samp$beta_6*val$Birth.Year*as.integer(val$TypeUsPres) + samp$beta_7*val$Birth.Year*as.integer(val$TypeChinaEmp) + samp$beta_8*val$Birth.Year*as.integer(val$TypeDalaiLama) + samp$beta_9*val$Birth.Year*as.integer(val$TypeJapanEmp)
    mu <- exp(beta_temp)
    t <- dweibull(samp$r, mu)
    ysims[i, k] <- t
  }
}

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

## Warning in dweibull(samp$r, mu): NaNs produced

r <- sample(nrow(non_censored_prediction), 1)
hist(ysims[r,])
abline(v=non_censored_prediction$Age.Event[r], col = "red")

print(mean(ysims<0))

## [1] NA

print(mean(non_censored_prediction$Age.Event <0))

## [1] 0

Lexis Diagram

#lg <- lexis_grid(1960, 2000, 30, 80, delta=5) 
#leaders_data$Initial.Age <- year(leaders_data$Event.Date) - leaders_data$Birth.Year
#subset <- filter(leaders_data, Censored == 0)
#subset <- subset[0:5,]
#ggplot() + geom_polygon(data = subset, mapping = aes(x = Birth.Year, y = Age.Event, group = Name)) + geom_point(data = subset, aes(x = Birth.Year, y = Birth.Year + Age.Event)) + geom_polygon(data = subset, mapping = aes(x = Birth.Year, y = ))

subset <- filter(leaders_data, Censored == 0)
ggplot(data = subset, aes(x = Birth.Year, Age.Event)) + geom_point(color="pink") + ggtitle("Distribution of Leaders' Ages") + labs(x = "Birth Year", y = "Age")